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ABSTRACT 

Context. Grains in circumstellar disks are believed to grow by mutual collisions and subsequent sticking due to surface forces. Results 
of many fields of research involving circumstellar disks, such as radiative transfer calculations, disk chemistry, magneto-hydrodynamic 
simulations largely depend on the unknown grain size distribution. 

Aims. As detailed calculations of grain growth and fragmentation are both numerically challenging and computationally expensive, 
we aim to find simple recipes and analytical solutions for the grain size distribution in circumstellar disks for a scenario in which 
grain growth is limited by fragmentation and radial drift can be neglected. 

Methods. We generalize previous analytical work on self-similar steady-state grain distributions. Numerical simulations are carried 
out to identify under which conditions the grain size distributions can be understood in terms of a combination of power-law distri- 
butions. A physically motivated fitting formula for grain size distributions is derived using our analytical predictions and numerical 
simulations. 

Results. We find good agreement between analytical results and numerical solutions of the Smoluchowski equation for simple shapes 
of the kernel function. The results for more complicated and realistic cases can be fitted with a physically motivated "black box" 
recipe presented in this paper. Our results show that the shape of the dust distribution is mostly dominated by the gas surface density 
(not the dust-to-gas ratio), the turbulence strength and the temperature and does not obey an MRN type distribution. 

Key words, accretion, accretion disks - protoplanetary disks - stars: pre-main-sequence, circumstellar matter - planets and satellites: 
formation 



1. Introduction 

Dust size distributions are a fundamental ingredient for many 
astrophysical models in the context of circumstellar disks 
and planet formation: whenever dust is present, it dominates 
the opacity of the disk, thereby influencing the temperature 
and consequently also the vertical structure of the disk (e.g., 
Dullemond & Dominik 2004). Small grains effectively sweep up 
electrons and therefore strongly affect the chemistry and the ion- 
ization fraction (also via grain surface reactions, see Vasyunin et 
al., in prep.) and thereby also the angular momentum transfer of 
the disk (e.g., Wardle & Ng 1999; Sano et al. 2000). 

Today, it is well established that the dust distributions in as- 
teroid belts and debris disks are governed by a so-called "col- 
lision cascade" (see Williams & Wetherill 1994): larger bodies 
, in a gas free environment exhibit high velocity collisions (> 
km s _1 ), far beyond their critical fragmentation threshold, which 
lead to cratering or even complete shattering of these objects. 
The resulting fragments in turn suffer the same fate, thus pro- 
ducing ever smaller grains down to sizes below about a few 
micrometers where Poynting-Robertson drag removes the dust 
particles (e.g., Wyatt et al. 1999). The grain number density dis- 
tribution in the case of such a fragmentation cascade has been 
derived by Dohnanyi (1969) and Williams & Wetherill (1994) 
and was found to follow a power-law number density distribu- 
tion n(m) oc rrT a with index a = ^ (which is equivalent to 
n(a) oc gT 3 5 ), with very weak dependence on the mechanical 
parameters of the fragmentation process. Tanakaetal. (1996), 
Makino et al. (1998) and Kobayashi & Tanaka (2010) showed 



that this result is exactly independent of the adopted collision 
model and that the resulting slope a is only determined by the 
mass-dependence of the collisional cross-section if the model of 
collisional outcome is self-similar (in the context of fluid dynam- 
ics, the same result was independently obtained by Hunt 1982 
and Pushkin & Aref 2002). The value of ^ agrees well with the 
size distributions of asteroids (see Dohnanyi 1969) and of grains 
in the interstellar medium (MRN distribution, see Mathis et al. 
1977; Pollack et al. 1985) and is thus widely applied, even at the 
gas-rich stage of circumstellar disks. 

However, in protoplanetary disks gas drag damps the mo- 
tions of particles. Very small particles are tied to the gas and, 
as a result, have a relative velocity low enough to make stick- 
ing feasible. The size distribution therefore deviates from the 
MRN power-law. Theoretical models of grain growth indicate 
that particles can grow to sizes much larger than a few /mi 
(see Nakagawa et al. 1981; Weidenschilling 1980, 1984, 1997; 
Dullemond & Dominik 2005; Tanaka et al. 2005; Brauer et al. 
2008a; Birnstiel et al. 2009). Indeed, the observational evidence 
suggests that growth up to cm-sizes is possible (e.g., Testi et al. 
2003; Natta et al. 2004; Rodmann et al. 2006; Ricci et al. 2010). 

However, as particles grow, they become more loosely cou- 
pled to the gas. This results in an increase in the relative ve- 
locity between the particles, a common feature of most sources 
of particles' relative velocity (i.e., turbulence, radial drift, and 
settling motions). Therefore, we expect that the assumption of 
perfect sticking will break down at some point and other colli- 
sional outcomes (bouncing, erosion, catastrophic disruption) be- 
come possible (see Blum & Wurm 2008). It is expected, then, 
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that at a certain point growth will cease for the largest particles 
in the distribution. Collisions involving these particles result in 
fragmentation, thus replenishing the small grains. On the other 
hand mutual collisions among small particles still result in co- 
agulation. As a result, a steady-state emerges. In this paper, this 
situation is referred to as a fragmentation-coagulation equilib- 
rium. 

The situation in protoplanetary disks differs, therefore, from 
that in debris disks. In the latter only fragmentation operates. 
The mass distribution still proceeds towards a steady state but, 
ultimately, mass is removed from the system due to, radiation 
pressure or Poynting-Robertson drag. 

In this paper, we consider the situation that the total mass 
budget in the system is conserved. For simplicity, we will ig- 
nore motions due to radial drift in this study. This mechanism 
effectively removes dust particles from the disk as well as pro- 
viding particles with a large relative motion. However, the de- 
rived presence of mm-size dust particles in protoplanetary disks 
is somewhat at odds with the usual (laminar) prescriptions for 
the radial drift rate (Weidenschilling 1977). Recently, it was 
shown that mm observations of protoplanetary disks can be ex- 
plained by steady-state size distributions if radial drift is inef- 
ficient (Birnstiel et al. 2010b). On the other hand, if radial drift 
would operate as effectively as the laminar theory predicts, then 
the observed populations of mm-sized particles at large disk radii 
cannot be sustained (Brauer et al. 2007). Possible solutions to 
reduce the drift rate include bumps in the radial pressure pro- 
file (see Kretke & Lin 2007; Brauer et al. 2008b; Cossins et al. 
2009) or zonal flows (see Johansen et al. 2009). 

In this work we analytically derive steady-state distributions 
of grains in the presence of both coagulation and fragmentation. 
The analytical predictions are compared to numerical simula- 
tions and applied to grain size distributions in turbulent circum- 
stellar disks. Both the theoretical and numerical results presented 
in this work are used to derive a fitting formula for steady-state 
grain size distributions in circumstellar disks. 

The paper is outlined as follows: in Sect. 2, we briefly 
summarize and then generalize previous results by Tanaka et al. 
(1996) and Makino et al. (1998). In Sect. 3, we test our theoret- 
ical predictions and their limitations by a state-of-the-art grain 
evolution code (see Birnstiel et al. 2010a). Grain size distribu- 
tions in circumstellar disks are discussed in Sect. 4. In Sect. 5, 
we present a fitting recipe for these distributions that can eas- 
ily be used in models where grain properties are important. Our 
findings are summarized in Sect. 6. 

2. Power-law solutions for a dust 

coagulation-fragmentation equilibrium 

In this section, we begin by summarizing some of the previ- 
ous work on analytical self-similar grain size distributions on 
which our subsequent analysis is based. We will then extend 
this to include both coagulation and fragmentation processes 
in a common framework. Under the assumption that the rel- 
evant quantities, i.e., the collisional probability between parti- 
cles, the distribution of fragments, and the size distribution, be- 
have like power-laws, we will solve for the size distribution in 
coagulation-fragmentation equilibrium. For simplicity, we con- 
sider a single monomer size only of mass niQ. This is therefore 
the smallest mass in the distribution. 

Another key assumption of our analytical model is that we 
assume the existence of a sharp threshold mass nif, above which 
collisions always result in fragmentation and below which colli- 



sion always result in coagulation. As explained above, the phys- 
ical motivation for this choice is that relative velocities increase 
with mass. This also means that in our theoretical model we ne- 
glect collision outcomes like bouncing or erosion (erosion is in- 
cluded in the simulations). Even though, small particles will in 
reality cause cratering/growth instead of complete fragmentation 
of the larger particle, we find that this assumption is often jus- 
tified because fragmenting similar-sized collisions prevent any 
growth beyond nif . Thus, collisions with much smaller particles 
(m « nif) do not have an important influence on the maximum 
size, however, they can significantly change the amount of small 
particles in the stationary distribution (see Sect. 5). We further 
assume a constant porosity of the particles, which relates the 
mass and size according to 

m = —Ps a . (1) 

where p s is the internal density of a dust aggregate. 

In our analysis, we will identify three different regimes, 
which are symbolically illustrated in Fig. 1: 

- Regime A represents a case where grains grow sequentially 
(i.e. hierarchically) by collisions with similar sized grains 
until they reach an upper size limit and fragment back to 
the smallest sizes. The emerging power-law slope of the size 
distribution depends only on the shape of the collisional ker- 
nel. 

- Regime B is similar to regime A, however in this case the 
fragmented mass is redistributed over a range of sizes and, 
thus, influencing the out-coming distribution. 

- In Regime C, the upper end of the distribution dominates 
grain growth at all sizes. Smaller particles are swept up by 
the upper end of the distribution and are replenished mostly 
by redistributed fragments of the largest particles. The result- 
ing distribution function depends strongly on how the frag- 
mented mass is distributed after a disruptive collision. 

For each of these regimes, we will derive the parameter ranges 
for which they apply and the slopes of the resulting grain size 
distribution. 

2.1. The growth cascade 

The fundamental quantity that governs the time-evolution of the 
dust size distribution is the collision kernel C mum . It is defined 
such that 

C mum2 ■ n{m x ) ■ n(m 2 ) dmidm 2 (2) 

gives the number of collisions per unit time per unit volume be- 
tween particles in mass interval [m\,m\ + dm\] and [ni2,ni2 + 
dni2], where n(m) is the number density distribution. Once spec- 
ified it determines the collisional evolution of the system. In the 
case that the number density n(m) does not depend on position, 
C mum2 is simply the product of the collision cross section and the 
relative velocity of two particles with the masses ;«i and ;«2. 

We use the same Ansatz as Tanaka et al. (1996), assuming 
that the collision kernel is given in the self-similar form 

C mi j n2 =m y 1 h\ — \. (3) 
\mi) 

Here, h is any function which depends only on the masses 
through the ratio of ;«2/mi. By definition, the kernel C mum2 has 
to be symmetric, therefore, Eq. 3 implies that h(ni\,ni2) is not 
symmetric (see Eq. 21). v is called the index of the kernel or the 
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Fig. 1: Illustration of the three different regimes for which ana- 
lytical solutions have been derived. Case A represents the growth 
cascade discussed in Sect. 2.1, case B the intermediate regime 
(see Sect. 2.2) and case C the fragmentation dominated regime 
which is discussed in Sect. 2.3. Particles are shattered once they 
reach the fragmentation barrier ;«f since collision velocities for 
particles > nif exceed the fragmentation threshold velocity. 



degree of homogeneity. As we will see in the following, v is one 
of the most important parameters determining the resulting size 
distribution. Different physical environments are represented by 
different values of v. Examples include the constant kernel (i.e., 
v — 0, mass independent), the geometrical kernel (i.e., v = 2/3, 
velocity independent) or the linear kernel (i.e., v = 1, as for 
grains suspended in turbulent gas). 

It is further assumed that the number density distribution of 
dust particles follows a power-law, 



n(m) — A ■ m 



(4) 



The time evolution of the mass distribution obeys the equa- 
tion of mass conservation, 



dmn(m) dF(m) 



at 



dm 



(5) 



where the flux F(m) does not represent a flux in a typical con- 
tinuous way since coagulation is non-local in mass space (each 
mass can interact with each other mass) but rather an inte- 
gration of all growth processes which produce a particle with 
mass greater than m out of a particle that was smaller than m 
(i.e. collisions of mi < m with any other mass m 2 such that 
m\ + mi > m). The flux in the case of pure coagulation was 
derived by Tanaka et al. (1996) 



F(m) 



dmi I 

mo J m- 



dm 2 mi C m . «_ n(mi) n(m 2 ), (6) 



where we changed the lower bound of the integration over mi to 
start from the mass of monomers mo (instead of 0) and the upper 
bound of the integral over m 2 to a finite upper end Wf (instead of 
infinity), compared to the definition used by Tanaka et al. (1996). 

Substituting the definitions of above and using the dimen- 
sionless variables Xi = mi/m, x 2 - m 2 /m, xq = mo/m, and 
Xf = ntf/m one obtains 



F(m) — rri 1 



-2a+3 



f djci f 

Jx„ Jl- 



d^2-^i 



"v,'7/|- 
x\) 



(7) 



where K approaches a constant value in the limit of m » mo and 
m <K »jf. 

Postulation of a steady state (i.e. setting the time derivative 
in Eq. 5 to zero), leads to the condition 



F(m) oc m v 2a+i = const., 
from which it follows that the slope of the distribution is 

v + 3 



(8) 



(9) 

This result was already derived for the case of fragmentation by 
Tanaka et al. (1996) and Dohnanyi (1969) and for the coagula- 
tion by Klett (1975) and Camacho (2001). The physical interpre- 
tation of this is a "reversed" fragmentation cascade: instead of a 
resupply of large particles which produces ever smaller bodies, 
this represents a constant resupply of monomers which produce 
ever larger grains (cf. case A in Fig. 1). 

2.2. Coagulation fragmentation equilibrium 

As Tanaka et al. (1996) and Makino et al. (1998) pointed out, 
the previous result is independent of the model of collisional 
outcomes as long as this model is self-similar (Eq. 3). However 
this is no longer the case if we consider both coagulation and 
fragmentation processes happening at the same time. 

We will now consider the case with a constant resupply of 
matter, due to the particles that fragment above mf. We assume 
these fragments obey a power-law mass distribution and are pro- 
duced at a rate 

h f (m) = N ■ m~t, (10) 

where £ reflects the shape of the fragment distribution and is 
a constant. 

If there is a constant flux of particles F(m) as defined above, 
then the flux of fragmenting particles (i.e., the flux produced by 
particles that are growing over the fragmentation threshold) is 
given by 

F(m f ) = K ■ m v f 2o+3 (11) 

where K is the integral defined in Eq. 7. 

The resulting (downward) flux of fragments Ff(m) can then 
be derived by inserting Eq. 10 into the equation of mass conser- 
vation (Eq. 5), 

dFf(m) 



dm 



— —m ■ rif. 



Integration from monomer size mo to m yields 



Ff(m) = -N ■ 



1 



2-f 



(m 2 ^ — m 



(12) 



(13) 



The normalization factor jV can be determined from the equilib- 
rium condition that the net flux vanishes, 



Ff(mf) — —F(mf), 



(14) 
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and was found to be 



N = (2 - £) K ■ 



^v-2a+3 



2-f 



(15) 



In Eq. 13, we can distinguish two cases: 



If £ > 2, the contribution of mo dominates the term in brack- 
ets. This means that most of the fragment mass is redis- 
tributed to monomer sizes and the situation is the same as in 
the pure coagulation case (cf. Case A in Fig. 1). The steady- 
state condition F(m) + Ff(m) = (i.e., the net flux is zero) 
yields that 



K ■ 



-N ■ 



1 



2-f 



2-f 



(16) 



is constant, which leads to the same result as Eq. 9. 
Intuitively, this is clear since the majority of the redistributed 
mass ends up at m ~ niQ. 

If £ < 2, the m-dependence dominates the term in brackets 
in Eq. 13 and postulation of a steady-state, 



K ■ m v - 2a+3 a N ■ 



1 



leads to an exponent 



2-? 
v + £+ 1 



(17) 



(18) 



less than Eq. 9. In this case, the slope of the fragment dis- 
tribution matters. This scenario is represented as case B in 
Fig. 1. 



2.3. Fragment dominated regime 

The result obtained in the previous section may seem to be quite 
general. However, it does not hold for low ^-values as we will 
show in the following. 

In our case, the integrals do not diverge due to the finite in- 
tegration bounds. However Makino et al. (1998) used and oo 
as lower and upper bounds for the integration and thus needed 
to investigate the convergence of the integral. They derived the 
following conditions for convergence: 



v — y — a + I < 
y - a + 2 > 0. 



(19) 
(20) 



where y gives the dependence of mijm\ in the /z-function of 
Eq. 21 (see Makino et al. 1998): 



mi 



h ■ - 



ife) 



for 2& « 1 



v-y 



(21) 



for a » 1 



The first condition (Eq. 19) considers the divergence towards 
the upper masses, whereas the second condition pertains the 
lower masses. We assume that Eq. 20 is satisfied and consider 
the case of decreasing £ for a given by Eq. 18 which results in a 
steeper size distribution (where the mass is concentrated close to 
the upper end of the distribution). We see that for £ < 1 + v - 2y, 
Eq. 19 is no longer fulfilled. The behavior of the flux integral 
changes qualitatively: growth is no longer hierarchical, but it be- 
comes dominated by contributions by the upper end of the inte- 
gration bounds. 



Physically, this means that the total number of collisions of 
any grain is determined by the largest grains in the upper end 
of the distribution. Hence, smaller sized particles are predomi- 
nantly refilled by fragmentation events of larger bodies instead 
of coagulation events from smaller bodies and they are predom- 
inantly removed by coagulation events with big particles (near 
the threshold ntf) instead of similar-sized particles. This corre- 
sponds to Case C in Fig. 1. 

To determine the resulting power-law distribution, we again 
focus on Eq. 6. The double integral of the flux F(m) can now be 
split into three separate integrals, 



■in, 



rf r> 

F(m) — I dm\ I d.m2m\C mumi n{m\)n{m2) 

Jmo Jm—mi 

Am\ I dm 2 m l C mum2 n(m l )n(m 2 ) (22) 

J m—mi 

+ I d/)?i I &m 1 m\C mumi n(m\)n{mi) 



according to whether ni2 is larger or smaller than m\. 

It can be derived (see Appendix A) that if the condition above 
(Eq. 19) is violated and if »jf » m, then the first and the third 
integral in Eq. 22 dominate the flux due to the integration until 
ntf. In this cases, the flux F(m) is proportional to m 2+y ~ a . 

A stationary state in the presence of fragmentation (Eq. 13) 
is reached if the fluxes cancel out, which leads to the condition 



a — ^ + y. 



(23) 



This is the sweep-up regime where small particles are cleaned 
out by big ones (cf. Case C in Fig. 1). 

2.4. Summary of the regimes 

Summarizing these findings, we find that the resulting distribu- 
tion is described by three scenarios (depicted in Fig. 1), depend- 
ing on the slope of the fragment distribution: 



Case A (growth cascade): 



£>2 



a - 



Case B (intermediate regime): v — 2y + 1 < ^ < 2 a 



v+3 



v+f+1 



(24) 



Case C (fragment dominated): £ < v — 2y + I a — £, + y 

3. Simulation results interpreted 

In this section, we will test the analytical predictions of the previ- 
ous section by a coagulation/fragmentation code (Birnstiel et al. 
2010a, see also Brauer et al. 2008a). The code solves for the time 
evolution of the grain size distribution using an implicit integra- 
tion scheme. This enables us to find the steady-state distribu- 
tion by using large time steps. In this way, the time evolution is 
not resolved, but the steady-state distribution is reliably and very 
quickly derived. 

We start out with the simplest case of a constant kernel and 
then - step by step - approach a more realistic scenario (in the 
context of a protoplanetary disk). In Sect. 4, we will then con- 
sider a kernel taking into account relative velocities of Brownian 
motion and turbulent velocities and also a fragmentation proba- 
bility which depends on the masses and the relative velocities of 
the colliding particles. 

The following results are only steady-state solutions, 
whether or not this state is reached depends on several condi- 
tions. Firstly, particles need to fragment. If there is no upper 
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boundary for growth, it will proceed unlimited and a steady state 
will never be reached. Secondly, radial motion needs to be slow 
enough to allow for a steady state. If quantities like the surface 
density or the temperature vary smoothly between neighbouring 
regions in the disk, the steady state solutions will also be simi- 
lar and radial transport will happen on top of a steady-state grain 
distribution. However if radial drift is acting strongly on particles 
of the distribution, the steady state will not be reached. Thirdly, a 
distribution of initially sub-//m sized grains will need some time 
to get into an equilibrium state. This time scale can be as small 
as < 1000 years inside of a few AU, while it can be of the order 
of a million years at 100 AU. We provide a rough estimate of 
this time scale in Appendix C. 



of £ = 0.5. For larger ^-values, the slope of the mass distri- 
bution flattens. In all cases, a bump develops towards the up- 
per end of the distributions. The reason for this "pile-up" is the 
following: grains typically grow mostly through collisions with 
similar-sized or larger particles. Since the distribution is trun- 
cated at the upper end (defined as a max , see also Eq. 43), particles 
close to the upper end lack larger collision partners, the growth 
rate at these sizes would decrease if the distribution keeps its 
power-law nature. This, in turn means that the flux could not be 
constant below a max . To keep a steady state, the number of par- 
ticles at that point has to increase in order to replace the missing 
collision partners at larger sizes. 



3.1. Constant kernel 

In the following section, we consider the case of a constant ker- 
nel and will include fragmentation above particle sizes of 1 mm 
because this represents an instructive test case. 




2.2 
2 
1.8 
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simulation result 
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Fig. 3: Exponent of grain size distributions for a constant kernel 
(i.e. v = 0) and different distributions of fragments (solid line) 
and the analytic solution for the growth cascade (Case A, dot- 
ted line), the intermediate regime (Case B, dashed line), and the 
fragment dominated regime (Case C, dash-dotted line). 



Fig. 2: Grain size distributions 1 for a constant kernel (i.e., v = 0) 
and different distributions of fragments. The peak towards the 
upper end of the distribution is due to the fragmentation barrier 
(explanation in the text). The slope of the mass distribution cor- 
responds to 6 - 3a. 



We iteratively solve for the steady-state size distribution be- 
tween coagulation and fragmentation. The outcome of these sim- 
ulations for a constant kernel (i.e., v = 0) are power-law distri- 
butions where the slope of the distribution depends on the frag- 
mentation law (the slope £ see Eq. 10). Figure 2 shows the cor- 
responding size distributions for some of the different fragment 
distributions: the steepest distribution corresponds to the case 

1 It should be noted that in this paper we will plot the distributions 
typically in terms of n(a) • m • a which is proportional to the distribution 
of mass. The advantage of plotting it this way instead of plotting n(a) 
is the following: when n(m) follows a power-law m' a , then the grain 
size distribution n(a) describes a power-law with exponent 2 - 3a and 
the mass distribution exponent is 6 - 3a. For typical values of a, this 
is less steep and differences between a predicted and a real distribution 
are more prominent. 



Figure 3 shows how the slope of the resulting distribution 
depends on the fragmentation slope f , where the three previously 
discussed regimes can be identified: 

- Case A, growth cascade: as predicted, this scenario holds for 
values of £ > 2 where most of the fragmenting mass is redis- 
tributed to fragments. This case corresponds to a "reversed" 
collisional cascade (just the direction in mass space is re- 
versed as collisions mostly lead to growth instead of shatter- 
ing). 

- Case B, the intermediate regime: when the fragment mass is 
more or less equally distributed over all sizes, mass gain by 
redistribution of fragments and the mass loss due to coagu- 
lation have to cancel each other. 

- Case C: most of the fragmented mass remains at large sizes. 
Therefore the mass distribution is dominated by the largest 
particles. In other words, growth is not hierarchical anymore. 
In this test case y — 0, and from Eq. 19 it follows that the 
transition between the intermediate and fragment-dominated 
regime lies at £ = 1, which can be seen in Fig. 3. 

The measured slopes of the size distribution for m«m[ are in 
excellent agreement with the model outlined in Sect. 2. 
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3.2. Including settling effects 

The simplest addition to this model is settling: as grains become 
larger, they start to settle towards the mid-plane. However, tur- 
bulent mixing counteracts this systematic motion. The vertical 
distribution of dust in a settling-mixing equilibrium can (close 
to the mid-plane) be estimated by a Gaussian distribution with a 
size-dependent dust scale height Smaller particles are well 
enough coupled to the gas to have the same scale height as the 
gas, H g , while larger particles decouple and their scale height is 
decreasing with grain size, 



(25) 




— = J — , for St > a t 



(see e.g., Dubrulle et al. 1995; Schrapler & Henning 2004; 
Youdin & Lithwick 2007) where the Stokes number St is a di- 
mensionless quantity which describes the dynamic properties of 
a suspended particle. Very small particles have a small Stokes 
number and are therefore well coupled to the gas. Particles which 
have different properties (e.g., size or porosity) but the same St 
behave aerodynamically the same. In our prescription of turbu- 
lence, St is defined as the product of the particles stopping time 
T st and the orbital frequency Q^. We focus on the Epstein regime 
where the Stokes number can be approximated by 



St = Ll k ■ T st - — ^ 



(26) 



with Eg being the gas surface density and p s being the inter- 
nal density of the particles which relates mass and size via 
m — 47r/3p s a 3 , 

Settling starts to play a role as soon as the Stokes number 
becomes larger than the turbulence parameter a t which can be 
related to a certain size 



2a t X g 



(27) 



Eq. 27 only holds within about one gas pressure scale height 
because the dust vertical structure higher up in the disk devi- 
ates from the a Gaussian profile (see also Dubrulle et al. 1995; 
Schrapler & Henning 2004; Dullemond & Dominik 2004). 

The mass-dependent dust scale height causes the number 
density distribution n(m) to depend on the vertical height, z. In 
disk-like configurations, it is therefore customary to consider the 
column density, 



N(m) 



n 

CO 



(m, z) dz. 



(28) 



Similar to Eq. 2, we can write the vertically integrated number 
of collisions as 



C mum2 ■ N(m\) ■ N(m 2 ) dm x dm 2 , 



(29) 



which gives the total number of collisions that take place over 
the entire column of the disks. 

The dependence of the collisional probability on scale 
height, is now reflected in the modified kernel C mumi (see 
Birnstiel et al. 2010a, Appendix A for derivation): 



r 



c„ 



^2* (fl? + flf) 



(30) 



The point to realize here is that, due to the symmetry be- 
tween Eq. 2 and Eq. 29, the analysis in Sect. 2 holds also for 



disk-like configuration, if the kernel is now replaced by C mum . 
The resulting exponent a then concerns the column density de- 
pendence (N(m) oc m~ a ). 

If we consider the case that St > a t and substitute Eq. 25 and 
Eq. 26 into Eq. 30, we find that 



r 



r 



^2k(h 2 1+ hI) 



= C 



Hi | Uj ' H j 



H 2 \ 
H 2 



mi 



(31) 



has an index v = 1 /6 for grain sizes larger than a sett and v = 
otherwise (it should be noted that Hi and H 2 are the dust scale 
heights whereas h(m*i_lm\) represents the function defined in 
Eq. 3). 

The theory described in Sect. 2 is strictly speaking only valid 
for a constant v-index, but if this index is constant over a signifi- 
cant range of masses, then the local slope of the distribution will 
still adapt to this index. In the case of settling, we can therefore 
describe the distribution with two power-laws as can be seen in 
Fig. 4. 

The fact that the distribution will locally follow a power-law 
is an important requirement for being able to construct fitting 
formulas which reproduce the simulated grain size distributions. 
It allows us in some cases to explain the simulation outcomes 
with the local kernel index (although coagulation and fragmenta- 
tion are non-local processes in mass space, since each mass may 
interact with each other mass). A physically motivated recipe to 
fit the numerically derived distribution functions for the special 
case of £ = 1 1/6 is presented in Sect. 5. 




10 



10 10 10 

grain size [cm] 



10" 



10" 



Fig. 4: Simulation result (solid line) and a "two power-law fit" 
to the vertically integrated dust distribution for a constant kernel 
with settling included. The dotted, vertical line denotes the grain 
size above which grains are affected by settling. 
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3.3. Non-constant kernels 

We performed the same tests as in Sect. 3.1 also for non- 
constant kernels, i.e. kernels with v > 0. The measured slopes 
for the cases of (v = 5/6, y = 0) (corresponding to the sec- 
ond turbulent regime, see Sect. 4.1), (v = 1/3, y = 1/3), and 
(v = 1/6,7 = ~l/2) (i.e., Brownian motion, see Sect. 4.1) are 
shown in Fig. 5. Similar to Fig. 3, the distribution always fol- 
lows the minimum index of the three different regimes (Cases 
A, B, and C). It should be noted that for all indices ^ of the frag- 
mentation law, all processes - coagulation, fragmentation, and 
re-distribution of fragments - take place, but the relative impor- 
tance of them is what determines the resulting slope. 

In Fig. 5, it can be seen that Case B (defined in Sect. 2.2) 
vanishes for the kernel in the upper panel, while it is present for 
a large range of £ for the kernel in the middle panel. This can 
be explained by the definitions of the three regimes which were 
summarized in Eq. 24: with (v = 5/6, y — 0) (cf. upper panel in 
Fig. 5), Case B is confined between £ = 11/6 and 2. The grain 
size distribution, therefore, switches almost immediately from 
being growth dominated (Case A) to fragmentation dominated 
(Case C). In the case of a kernel with v = 1/3 and y = 1/3, this 
range spans from 2/3 to 2, as can be seen in the central panel of 
Fig. 5. 

The lower panel shows the distribution for a Brownian mo- 
tion kernel (i.e., v = 1/6 and y = -1/2). The grey shaded area 
highlights the range of £ values where our predictions do not ap- 
ply, that is, Eq. 20 is no longer fulfilled and thus, the resulting 
steady-state distributions are no longer power-law distributions. 



v = 5/6, 7=0 



0.5 - 




simulation result 
growth cascade (Case A) 
intermediate regime (Case B) 
fragment dominated (Case C) 



v = 1/3, 7 = 1/3 




v = 1/6, 7 = -1/2 



4. Grain size distributions in circumstellar disks 

In this section, we will leave the previous "clean" kernels and 
focus on the grain size distribution in circumstellar disks includ- 
ing relative velocities due to Brownian and turbulent motion of 
the particles, settling effects, and a fragmentation probability as 
function of particle mass and impact velocity. 

Combining all these effects makes it impossible to find 
simple analytical solutions in the case of a coagulation- 
fragmentation equilibrium. We will therefore use a coagulation 
fragmentation code to find the steady-state solutions and to show 
how the steady-state distributions in circumstellar disks depend 
on our input parameters. 

We will first discuss the model "ingredients", i.e. the rela- 
tive velocities and the prescription for fragmentation/sticking. 
Afterwards, we will define characteristic sizes at which the shape 
of the size distribution changes due to the underlying physics. In 
the last subsection, we will then show the simulation results and 
discuss the influence of different parameters. 

4.1. Relative velocities 

We will now include the effects of relative velocities due to 
Brownian motion, and due to turbulent mixing. The Brownian 
motion relative velocities are given by 



A«BM - 



8kbT (mi + 1112) 
nm\ ni2 




2.5 



£ where hf(m) oc m 



Fig. 5: Exponents of the grain size distributions for three differ- 
ent kernels as function of the fragment distribution index. The 
plot shows the simulation results (solid lines), the analytic solu- 
tion for the growth cascade (Case A, dotted lines), the interme- 
diate regime (Case B, dashed lines), and the fragment dominated 
regime (Case C, dash-dotted lines). The upper panel was calcu- 
lated with a (v = 5/6, y = 0)-kernel (i.e. turbulent relative veloc- 
ities, see Sect. 4.1), the middle panel with a (v = 1/3, y = 1/3)- 
kernel and the lower panel with a (v = 1/6, y = -l/2)-kernel 
(i.e. Brownian motion relative velocities). In the grey shaded 
area, Eq. 20 is not fulfilled, and the size distribution is not a 
power-law. 



(32) 



where k\> is the Boltzmann constant and T the mid-plane tem- 
perature of the disk. Ormel & Cuzzi (2007) have derived closed 
form expressions for particle collision velocities induced by tur- 
bulence. They also provided easy-to-use approximations for the 



different particle size regimes which we will use in the follow- 
ing. Small particles (i.e. stopping time of particle <k eddy cross- 
ing time) belong to the first regime of Ormel & Cuzzi (2007) 
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Table 1 : Kernel indices for the different regimes without settling. 



Regime 


V 


r 


upper end 


Brownian motion regime 


i 

6 


i 

2 


fl BT 


Turbulent regime I 


1 





a ,2 


Turbulent regime II 


5 
6 





"max 



with velocities proportional to 

Ami k [Sti - St 2 | . 



(33) 



The relative velocities of particles in the second turbulent regime 
of Ormel & Cuzzi (2007) are given by 



Amti 



(34) 



where St max is the larger of the particles Stokes numbers. 
Velocities in this regime show also a weak dependence on the 
ratio of the Stokes numbers which we will neglect in the follow- 
ing discussion. 

Together with the geometrical cross section <x geo = ji(a\ + 
a 2 ) 2 , it is straight-forward to estimate the indices of the ker- 
nel, v and y, as defined in Eq. 3 and 21 for all these regimes, 
without settling (assuming that only Brownian motion or turbu- 
lent motion dominates the relative velocities). The indices for 
these three sources of relative particle motion are summarized 
in Table 1 (for a derivation, see Appendix B). If settling is to be 
included, then the v index for particle sizes above a sett has to be 
increased by g (see Sect. 3.2) while y remains the same. The v 
and y indices for all three regimes can be found in Table 1. 



4.2. Fragmentation and watering 

We introduce fragmentation and cratering according to the 
recipe of Birnstiel et al. (2010a): whether a collision leads to 
sticking or to fragmentation/cratering is determined by the frag- 
mentation probability 



Pi = 







1 - 



if Am < Mf • 
if Au > Uf 



5u 



(35) 



^ else 



which means that all impacts with velocities above the critical 
break-up threshold Mf lead to fragmentation or cratering while 
impacts with velocities below Mf - 5u lead to sticking. The width 
of the linear transition region 5u is chosen to be 0.2 Mf, since 
laboratory experiments suggest that there is no sharp fragmen- 
tation velocity (Blum & Muench 1993; C. Guttler 2010, private 
communication). Our simulation results do not show a strong 
dependence on this value, which was tested by varying 5u by a 
factor of 2. 

The mass ratio of the two particles determines whether the 
impact completely fragments the larger body (masses within one 
order of magnitude) or if the smaller particle excavates mass 
from the larger body (masses differ by more than one order of 
magnitude). This distinction between an erosion and a shatter- 
ing regime follows the numerical studies of Paszun & Dominik 
(2009) and Ormel et al. (2009) and experimental studies of 
Guttler et al. (2010). These works do not precisely constrain the 
mass ratio which distinguishes between both regimes, however 
our simulation results do only weakly depend on it. 



In the case of fragmentation, the whole mass of both colli- 
sion partners is redistributed to all masses smaller than the larger 
body according to Eq. 10. 

In the case of cratering, it is assumed that the smaller particle 
(with mass m imp ) excavates its own mass from the larger body. 
The mass of the impactor as well as the excavated mass is then 
redistributed to masses smaller than the impactor mass according 
to Eq. 10. Thus, the total redistributed mass equals 



c 1, 

Jmn 



n(m) ■ m Am — 2 m 



imp 



(36) 



and the mass of the larger body is reduced by 2 m lmp . 

Most parameters such as the fragmentation velocity or the 
amount of excavated material during cratering are not yet well 
enough constrained (for the most recent experimental work, 
see Blum & Wurm 2008, Guttler et al. 2010, and references 
therein). Experiments suggest fragmentation velocities of a few 
m s _1 and fragment distributions with ^ between 1 and 2 
(Guttler et al. 2010 find £, values between 1 .07 and 1 .37 for Si0 2 
grains). Simulations of silicate grain growth around 1 AU show 
that also bouncing (i.e. collisions without sticking or fragmen- 
tation) can play an important role (see Weidling et al. 2009; 
Zsom et al. 2010). However, changes in material composition 
such as organic or ice mantles or the monomer size are expected 
to change this picture. As there is still a large parameter space to 
be explored, we continue with the rather simple recipe of stick- 
ing, fragmentation and cratering outlined above. 

4.3. Regime boundaries 

From Eq. 24, we can calculate the slope of the distribution in the 
different regimes if we assume that the slope of the distribution 
at a given grain size always follows from the kernel index v. To 
construct a whole distribution consisting of several power-laws 
for each regime, we need to know where each of the different 
relative velocity regime applies. 

It is important to note that relative velocities due to Brownian 
motion decrease with particle size whereas the relative veloci- 
ties induced by turbulent motion increase with particle size (up 
to St = 1). Therefore, Brownian motion dominates the relative 
velocities for small particles, while for larger particles, turbu- 
lence dominates. From numerical simulations, we found that at 
those sizes where the highest turbulent relative velocities (i.e. 
collisions with the smallest grains) start to exceed the smallest 
Brownian motion relative velocities (i.e. collisions with similar 
sized grains), the slope of the distribution starts to be determined 
by the turbulent kernel slope. By equating the approximate rela- 
tive velocity of Ormel & Cuzzi (2007) and Eq. 32, the according 
grain size can be estimated to be 



flBT 



8S g 



Re 



H m p ( An 
3n a t 



-Pi 



(37) 



where we approximate the Reynolds number (i.e., the ratio of 
turbulent viscosity v t = ac s H p over molecular viscosity) near 
the disk mid-plane by 



Re 



fft 2 g cr H , 
2 /i nip 



(38) 



Here, cr H , is the cross section of molecular hydrogen (taken to 
be 2 x 10~ 15 cm 2 ) and fi = 2.3 is the mean molecular weight in 
proton masses m p . 
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grain size [cm] grain size [cm] 

Fig. 6: Fiducial model and variations of the most important parameters: fragmentation power-law index threshold fragmentation 
velocity Uf, turbulence parameters a t , surface density E g , particle internal density p s , and mid-plane temperature T. The shape of the 
vertically integrated size distributions does not depend on the stellar mass or the distance to the star (only via the radial dependence 
of the parameters above). 



Turbulent relative velocities strongly increase for grains with 
a stopping time that is larger or about the turn-over time of the 
smallest eddies. More specifically, the Stokes number of the par- 
ticles at this change in the relative velocity is 

1 t„ 1 1 
Sti 2 = — = -Re", (39) 

where fL = 1/Ok and t n — ?l • Re" are the turn-over times of the 
largest and the smallest eddies, respectively, and Qk is the Kepler 



frequency. Ormel & Cuzzi (2007) approximated the factor y a to 
be about 1.6. The corresponding grain size in the Epstein regime 
is therefore given by 

1 2S 

a n = ^-Re^, (40) 

As mentioned above, the Brownian motion relative velocities of 
small grains decrease with their size. For larger sizes, the rela- 
tive velocities due to turbulent motion are gaining importance, 
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which are increasing with size until a Stokes number of unity. 
For typical values of the sound speed 

e 8 = M (41) 
\Jfim p 

and the turbulence parameter a t , the largest turbulent relative 
velocity AM max w Va^c s exceed the critical collision velocity 
of the grains (which is of the order of a few m s _1 ) and therefore 
leads to fragmentation of the dust particles. In the case of very 
quiescent environments and/or larger critical collision velocities, 
particles do not experience this fragmentation barrier and can 
continue to grow. Hence, a steady state is never reached. The 
work presented here focuses on the former case where 

AM max > Uf, (42) 

and grain growth is always limited by fragmentation. 

As relative turbulent velocities are (in our case) increas- 
ing with grain size, we can relate the maximum turbulent rel- 
ative velocity and the critical collision velocity to derive the ap- 
proximate maximum grain size which particles can reach (see 
Birnstiel et al. 2009) 

^max ~ ' ~~ T ■ (43) 

4.4. Resulting steady-state distributions 

The parameter space is too large to even nearly discuss all pos- 
sible outcomes of steady state grain size distributions. We will 
therefore focus on a few examples and rather explain the basic 
features and the most general results only. For this purpose, we 
will adopt a fiducial model and consider the influences of several 
parameters on the resulting grain size distribution: Wf, a t , S g , 
p s , and T (see Fig. 6). 

The fiducial model (see the solid black line in Fig. 6) shows 
the following features: steep increase from the smaller sizes un- 
til a few tenth of a micrometer. This relates to the regime dom- 
inated by Brownian motion relative velocities. The upper end 
of this regime can be approximated by Eq. 37. The flatter part 
of the distribution is caused by a different kernel index v in the 
parts of the distribution which are dominated by turbulent rel- 
ative velocities. The dip at about 60 /urn (cf. Eq. 40) is due to 
the jump in relative velocities as the stopping time of particles 
above this size exceeds the shortest eddy turn-over time (see 
Ormel & Cuzzi 2007). 

The upper end of the distribution is approximately at a max . 
The increased slope of the distribution and the bump close to the 
upper end are caused by two processes. Firstly, a boundary ef- 
fect: grains mostly grow by collisions with similar or larger sized 
particles. Grains near the upper end of the distribution lack larger 
sized collision partners and therefore the number density needs 
to increase in order to keep the flux constant with mass (i.e. in 
order to keep a steady-state). Secondly, the bump is caused by 
cratering: impacts of small grains onto the largest grains do not 
cause growth or complete destruction of the larger bodies, in- 
stead they erode them. Growth of these larger bodies is therefore 
slowed down and, similar to the former case, the mass distribu- 
tion needs to increase in order to fulfill the steady-state criterion 
("pile-up effect"). 

The upper left panel in Fig. 6 shows the influence of the 
distribution of fragments after a collision event: larger values 
of £ mean that more of the fragmented mass is redistributed 



to smaller sizes. Consequently, the mass distribution at smaller 
sizes increases relative to the values of smaller £ values. 

The strong influence of the fragmentation threshold velocity 
Uf can be seen in the upper right panel in Fig. 6: according to 
Eq. 43, an order of magnitude higher Uf leads to a 100 times 
larger maximum grain size. 

The grain size distributions for different levels of turbulence 
are shown in the middle left panel of Fig. 6. The effects are two- 
fold: firstly, an increased a t leads to increased turbulent relative 
velocities, thus, moving the fragmentation barrier a max to smaller 
sizes (cf. Eq. 43). Secondly, a larger a t shifts the transition within 
the turbulent regime, 012, to smaller sizes. Consequently, the sec- 
ond turbulent regime gains importance as a t is increased since its 
upper end lower boundary extend ever further. 




grain size [cm] 



Fig. 7: Comparison of the turbulent relative velocities of 
Ormel & Cuzzi (2007) to the fitting formula used in the recipe 
(see Eqs. 48 and 45) for grain size distributions. The largest er- 
ror in the resulting upper grain size op derived from the fitting 
formula is about 25%. 



The middle right panel in Fig. 6 displays the influence of an 
decreased gas surface density S g (assuming a fixed dust-to-gas 
ratio). It can be seen that not only the total mass is decreased due 
to the fixed dust-to-gas ratio but also the upper size of the dis- 
tribution a max decreases. This is due to the coupling of the dust 
to the gas: with larger gas surface density, the dust is better cou- 
pled to the gas. This is described by a decreased Stokes number 
(see Eq. 26) which, in turn, leads to smaller relative velocities 
and hence a larger a max . Interestingly, the shape of the grain size 
distribution does not depend on the total dust mass, but on the 
total gas mass, as long as gas is dynamically dominating (i.e., 
Z g » Sd). If more dust were to be present, grains would collide 
more often, thus, a steady-state would be reached faster and the 
new size distribution would be a scaled-up version of the former 
one. However the velocities at which grains collide are deter- 
mined by the properties of the underlying gas disk. In this way, 
the dust grain size distribution is not only a measure of the dust 
properties, but also a measure of the gas disk physics like the gas 
density and the amount of turbulence. 

The shape of the size distribution for different grain volume 
densities p s does not change significantly. However most regime 
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Table 2: Parameter values for which the recipe presented in 
Sect. 5 has been compared to simulation results: a A is the tur- 
bulence parameter, T is the mid-plane temperature, £ e is the gas 
surface density, and Uf is the critical collision velocity. 



parameter 


unit 






values 


a, 




io- 4 


io- 3 


IO" 2 - 


T 




10 


100 


500 10 3 - 


Eg 


[g cm" 2 ] 


0.1 


1 


10 100 10 3 10 4 


Uf 


[m s" 1 ] 


1 


3 


10 



boundaries (a se tt, <Zi2, and a max ) are inversely proportional to p s 
(because of the coupling to the gas, described by the Stokes num- 
ber). A decrease (increase) in p s therefore shifts the whole dis- 
tribution to larger (smaller) sizes as can be seen in the lower left 
panel in Fig. 6. 

The upper end of the distribution, a max , is inversely propor- 
tional to the mid-plane temperature T (as in the case of the tur- 
bulence parameter) whereas the transition between the different 
turbulent regimes an does not. Therefore, increasing the tem- 
perature decreases a max in the same way as decreasing E g does. 
However an is not influenced by temperature changes, there- 
fore the shape of the size distribution changes in a different way 
than in the case of changing Z g as can be seen by comparing the 
middle right and the lower right panels of Fig. 6. 



5. Fitting formula for steady-state distributions 

In this section, we will describe a simple recipe which allows us 
to construct vertically integrated grain size distributions which 
fit reasonably well to the simulation results presented in the pre- 
vious section. 

The recipe does not directly depend on the radial distance 
to the star or on the stellar mass. A radial dependence only en- 
ters via radial changes of the input parameters listed in Table 4. 
This recipe has been tested for a large grid of parameter values 2 
(shown in Table 2), however there are some restrictions. 

5.1. Limitations 

These fits strictly apply only for the case of £ = 11/6. In this 
case, the slopes of the distribution agree well with the predic- 
tions of the intermediate regime (Case B, defined in Sec. 2.2). 
For smaller values of £ the slopes do not strictly follow the an- 
alytical predictions. This is due to the fact that we include cra- 
tering, which is not covered by our theory. Erosion is therefore 
an important mode of fragmentation: it dominates over complete 
disruption through the high number of small particles (see also 
Kobayashi & Tanaka 2010) and it is able to redistribute signifi- 
cant amounts of mass to the smallest particle sizes. 

One important restriction for this recipe is the upper size of 
the particles a max : it needs to obey the condition of Eq. 42, since 
otherwise, particles will not experience the fragmenting high- 
velocity impacts and a steady state will never be reached since 
particles can grow unhindered over the meter-size barrier. 

There are also restrictions to a very small a max : if a max is 
close to or even smaller than 012, then the fit will not represent 
the true simulation outcome very well. In this case, the upper 
end of the distribution, and in more extreme cases the whole 
distribution will look much different. Thus, the sizes should obey 



the condition 

5pm < an < a max , 
where each inequality should be within a factor of a few. 



(44) 



5.2. Recipe 

The following recipe calculates the vertically integrated mass 
distribution of dust grains in a turbulent circumstellar disk within 
the the above mentioned limitations. The recipe should be ap- 
plied on a logarithmic grain size grid a, with a lower size limit 
of 0.025 pm and a fine enough size resolution (aj+i/flj < 1.12). 
For convenience, all variables are summarized in Table 4. The 
steps to be performed are as follows: 

1. Calculate the grain sizes which represent the regime bound- 
aries asj, a\2, flsett which are given by Eqs. 37, 40, and 27. 

2. Calculate the turbulent relative velocities for each grain size. 
For this, we approximate the equations which are given by 
Ormel & Cuzzi (2007). Collision velocities with monomers 
are given by 



Re J • (St, - St ) 



Am" 



*gas 



for at < a 12 



(1 - e) ■ Re* • (St,- - Sto) for an < a; < 5 an 
+ e ■ V3~St~ 



V3~st; 



for a,- > 5 fli 



(45) 

where Re is the Reynolds number (see Eq. 38), and St, 
and Sto are the Stokes numbers of a, and monomers (ao = 
0.025/mi), respectively (cf. Eq. 26). w gas is given by 



■■gas 



and the interpolation parameter e is defined as 



e = 



a; - a 12 

4fli 2 



(46) 



(47) 



and collisions with similar sized bodies are approximated as 
for a, < an 



Am 



VI 



(48) 



• Am™" for a* > a n 



See also www.mpia.de/distribution-fits 



A comparison between these approximations and the formu- 
las of Ormel & Cuzzi (2007) is shown in Fig. 7. 
Using Eqs. 45 and 48 for the relative velocities, and the tran- 
sition width 611 = 0.2Mf, find the grain sizes which corre- 
spond to the following conditions: 

- a^: particles above this size experience impacts with 
monomers with velocities of Am™ 011 > Uf - 6u (i.e., cra- 
tering starts to become important). 

- ap: particles above this size experience impacts with 
equal sized grains with velocities of AM^ q > Uf - 6u (i.e., 
fragmentation becomes important). 

- or: particles above this size experience impacts with sim- 
ilar sized grains with velocities of Au^ > Uf (i.e. every 
impact causes fragmentation/cratering). 
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Fig. 8: Step-by-step construction of the fit distribution for the following parameters: S g = 20 g cm 2 , Zd = 0.2 g cm 2 , a t = 1x10 4 , 
Uf - 1 m s _1 , £ = 1.833, T — 50 K, p s — 1.6 g cm -3 . 77je upper left panel shows the distribution after step 5: each interval obeys 
a different power-law, the distribution is continuous apart from a jump at an- The upper right panel displays the fit including an 
increase at the upper end, according to step 6. The bump caused by cratering (cf. Eq. 53) is shown in the bottom left panel while 
the bottom right panel compares the final fit distribution (solid curve) to the simulation result (dashed curve). The vertical lines 
correspond to the regime boundaries: «bt (solid), a se tt (dash-dotted), an (dashed) and a P (dotted). 



4. Calculate the factor J according to the recipe 



V = c s 



jump Eg \ 4 3 a t 



2.5" 9 + 



1.1 9 + 1+2 V3 



V 



Uf/ j 



9\ 



■1\ 



(49) 



(50) 



where p. = 2.3, cr H , = 2x 10" 15 cm 2 andj a = 1.6. 
The power-law indices of the mass distribution d>,- for each in- 
terval between the regime boundaries («bt, an, <hea, a ?) are 
calculated according to the intermediate regime (cf. Eq. 18). 
The slopes have to be chosen according to the kernel regime 
(Brownian motion, turbulent regime 1 or 2) and according 
to whether the regime is influenced by settling or not (i.e., 
if a is larger or smaller than a se tt)- The resulting slopes of 
the mass distribution are given in Table 3. The first version 
of the fit ficii) (where a, denotes the numerical grid point of 



the particles size array) is now constructed by using power- 
laws (oc a.') in between each of the regimes, up to op. The fit 
should be continuous at all regime boundaries except a drop 
of 1/J at the transition at an. An example of this first version 
is shown in the top left panel of Fig. 8. 
Mimic the cut-off effects which cause an increase in the dis- 
tribution function close to the upper end: linearly increase 
the distribution function for all sizes a mc < a < ap: 



f(fld->f(ad- 2 



-a P 



(51) 



with 



a inc = 0.3 a P . (52) 

The resulting fit after this step is shown in the upper right 
panel of Fig. 8. 
7. The bump due to cratering is mimicked by a Gaussian, 



b(ad = 2 • f(a h ) ■ exp 



(a,- - a P ) 



(T- 



(53) 
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Table 3: Power-law exponents of the distribution n(m) ■ m ■ a. 
The slopes were calculated using the formula for a coagula- 
tion/fragmentation equilibrium (Eq. 18). Within each regime of 
relative velocities, it has to be differentiated whether grains are 
influenced by settling or not. 



Regime 


Si 




a, 


> "sett 


Brownian motion regime 


3 
2 




5 
4 


Turbulent regime I 


1 

4 







Turbulent regime II 


1 

2 




i 

4 



where cr is defined as 

min(|a R - a P \ , \a L - a P [) 



cr = 



VM2) 

but should be limited to be at least 
cr > 0.1 flp. 



(54) 



(55) 



8. The fit T(a\) is now constructed by using the maximum of 
f(flt) and b{ai) in the following way: 



f(ad 



if a -i < «l 



max (/(a;), b(ai)) if < at < ap 
b(cii) if ap < di < «r 



(56) 







else 



Finally, the fit needs to be normalized to the dust surface 
density at the given radius. T is a (yet un-normalized) verti- 
cally integrated mass distribution (shown in the bottom right 
panel of Fig. 8). To translate this mass distribution to a ver- 
tically integrated number density distribution N(a), we need 
to normalize it as 



N(a) 



T(a) 



f °° T(o) dlna m-a 

Jan 



(57) 



6. Conclusions 

In this work, we generalize the analytical findings of previ- 
ous works to the case of grain size distributions in a coagula- 
tion/fragmentation equilibrium. Under the assumption that all 
grains above a certain size, a max , fragment into a power-law dis- 
tribution rtf(m) oc m~£, we derived analytical steady-state solu- 
tions for self-similar kernels and determined three different cases 
(see 2.4). Cratering is not covered by our theory. However our 
simulations that include cratering agree with the theoretical pre- 
dictions for a fragmentation law with £ = 11/6. 

Results show that dust size distributions in circumstel- 
lar disks do not necessarily follow the often adopted MRN 
power-law distribution of n(a) oc a~ 35 (see Mathis et al. 1977; 
Dohnanyi 1969; Tanaka et al. 1996; Makino et al. 1998; Garaud 
2007) when both coagulation and fragmentation events operate. 
We performed detailed simulations of grain growth and frag- 
mentation to test the analytical predictions and found very good 
agreement between the theory and the simulation results. 



We applied the theory to the gaseous environments of cir- 
cumstellar disks. Unlike the models of Garaud (2007), the upper 
end of the size distribution is typically not limited by the growth 
time scale but by fragmentation because relative velocities in- 
crease with grain size and reach values large enough to frag- 
ment grains. The shape of the dust distribution is determined 
by the gaseous environment (e.g., gas surface density, level of 
turbulence, temperature and others) since the gas is dynamically 
dominant as long as the gas surface density significantly exceeds 
the dust surface density. The total dust mass merely provides the 
normalization of the distribution and the time scale in which an 
equilibrium is reached. The results presented in this work show 
that the physics of growth and fragmentation directly link the 
upper and the lower end of the dust distribution in circumstellar 
disks. 

A ready-to-use recipe for deriving vertically integrated dust 
size distributions in circumstellar disks for a fixed value of 
£ = 11/6 is presented in Sect. 5. Although the collision kernel 
in circumstellar disks is complicated, we found good agreement 
with our fitting recipe for a fragment distribution with £ = 11/6. 
The recipe can readily be used for further modeling such as disk 
chemistry or radiative transfer calculations. 
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Table 4: Definition of the variables used in this paper. The variables grouped as "input variables" are the parameters of the fitting 
recipe in Sect. 5. "Other variables" summarizes the definitions of all other variables used in this recipe. 





variable 


definition 


unit 


03 




Turbulence strength parameter 


- 


1 




Fragmentation threshold velocity 


cm s~' 






Gas surface density 


g cm" 2 


> 

-4—1 


Id 


Dust surface density 


g CUT 2 


3 
& 


f 


Power-law index of the mass distribution of fragments, see Eq. 10 


_ 


d 


T 

Ps 


Mid-plane temperature 

Volume density of a dust particle 


K 

g cm" 3 




a 


grain size, a = (3 m/(4^p s )) 1/3 


cm 




a 


Slope of the number density distribution n(m) oc nT a , see Eq. 4 


_ 




fl B T 


Eq. 37 


cm 




fll2 


Eq. 40 


cm 




^max 


Eq. 43 


cm 




fl L 


Left boundary of the bump function b(a), see Sect. 5, paragraph 3 


cm 




flp 


Peak size of the bump function b(a), see Sect. 5, paragraph 3 


cm 




«R 


Right boundary of the bump function b(a), see Sect. 5, paragraph 3 


cm 




^sett 


Eq. 27 


cm 




"inc 


Eq. 52 


cm 




tXa) 


Bump function, see Eq. 53 






Cm i, mi 


Collision kernel, see Eq. 3 


3 -1 

cm s 


<X 
(0 


r 

y ~ , m\ ,m-> 


Collision kernel including settling effects, see Eq. 3 


2 -1 

cm s 


1 


c s 


sound speed, see Eq. 41 


cm s" 1 


'§ 


Si 


Slopes of the fit-function, see Table 3 




Other v 


Su 


Width of the transition between sticking and fragmentation, taken to be 0.2 Uf 


cm s" 1 


e 


Interpolation parameter, see Eq. 47 




7 
J 
K 


Power-law index of the function h(m 2 /mi) for large ratios of m 2 /in as defined in Eq. 21 
Empirical parametrization of the discontinuity, see Eq. 50 
Integral defined in Eq. 7 


_ 
_ 




k b 


Boltzmann constant 


erg K" 1 




m 


mass of the particle, m = 4tt/3 p s a 3 


g 




m 


Monomer mass 


g 




m i 


Mass above which particles fragment 


g 




1* 


Mean molecular weight in proton masses, taken to be 2.3 






m v 


Proton mass 


g 




n(m) 


Number density distribution 


g cm -3 




N(m) 


Vertically integrated number density distribution 


g cm" 2 




V 


Degree of homogeneity of the kernel as defined in Eq. 3 


_ 




Re 

P s 


Reynolds number, see Eq. 38 
density of a dust grain 


- 

g cm" 3 




St 


Mid-plane Stokes number in the Epstein regime, see Eq. 26 






0"H 2 


Cross-section of molecular hydrogen 


cm- 






mean square turbulent gas velocity, see Eq. 46 


cm s" 1 






Relative velocities between monomers and grains of size a,, see Eq. 45 


cm s" 1 




Au eq 

I 


Relative velocities between grains of size a,, see Eq. 48 


cm s" 1 




V 


Parameter used in the calculation of J, see Eq. 49 


cm s" 1 






equals 1.6, parameter from Ormel & Cuzzi (2007) 
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Appendix A: Derivation of the fragment dominated 
size distribution 

In this section, we will derive the slope of the size distribution 
which is dominated by the largest particles. Since in our sce- 
nario, the integrals are confined between the monomer mass mo 
and the largest particles at the fragmentation barrier »jf, the in- 
tegrals do not diverge as in the scenario of Tanaka et al. (1996) 
and Makino et al. (1998), who consider integration bounds of 
zero and infinity. However, if the first condition of Makino et al. 
(1998) is not fulfilled, then the mass flux (cf. Eq. 6) is dominated 
by the upper bound of the integral. This is the case which we will 
consider in the following. 

As noted in Sect. 2.3, the flux integral (Eq. 6) can be split into 
three separate integrals, in order to distinguish cases of m% > mi 
or m.2 < mi, 

F(m) = h + I 2 + h, (A.l) 

where 7i, I2, and 73 correspond (from left to right) to the three 
integrals defined in Eq. 22. We will now evaluate these integrals 
in the limits of mo <sc m <sc ntf, using the limiting behavior of 
h(milmi) as given in Eq. 21. 



The term in square brackets can be split into a sum of integrals, 
the first of which is straight-forward to evaluate as 



h 
D 



1 



(i) 



v - 2a + 3 



r, 1 + 

■ I axi x l 



v — y — a 



- a + 1 



(A.8) 



while the second can be identified as a sum of a Beta function 
B(a, b) and an incomplete Beta function Bi(a, b), 



h 
D 



1 



a+b-l 



a + b-l 
The conditions from above, 



(B(a,b)- B,(a,b)). (A.9) 



a-2+v-y-a > 
b = 2 + y - a > 



(A.10) 
(A.ll) 



assure that the Beta functions, and thus I2, are real. The numer- 
ical value of I2/D for a range of values of a and b are shown in 
Fig. A.l. 



A. 1. First integral 

Carrying out the integration over m2 in Ii leads to 

l2 . nm/2 



v-y-a+1 J m 

f v—y—a+ 1 



drai m 



y-a+ 1 



— (m — mi) 



v— y— a+\ 



]■ 



(A.2) 



from which we can derive Eq. 19, the first convergence criterion 
of Makino etal. (1998). 

We consider the case where this condition does not hold, i.e., 



v — y — a + I >0. 



(A.3) 



Then, the ntf term in brackets dominates over the other term. 
Thus, the term in brackets is constant and carrying out the inte- 
gration yields 



a 2 1 v—y-a+ 1 

A L ■ ho ■ m f " 
(v - y — a + 1) • (y — a + 2) 



(?) 



y-a+2 



y-a+2 



(A.4) 



Now, if the second condition, Eq. 20 holds, using mo <k m, we 
derive 



a 2 j v-y-a+l 

A ■ ho ■ m f ' 

II — 

(v — y — a + 1) • (y — a + 2) 
A.2. Second integral 



(?) 



y-a+2 



(A.5) 



We rewrite I2 using the dimensionless variables xi = mi/m and 
X2 - m2/m which yields 



12 — A ■ ho ■ m 



J I Jl-x, 



1 l+v-y-a y-a 

ax2 x 1 • x 7 . 



(A.6) 



By integrating over x%, we derive 
A 2 -h - m 3+v - 2a 



h 



y — a + 



t-v— 2a r*l 



1 +v — y — ct \ y— a+l 

axix l ■ bCj 



:=D 



(A.7) 




5 b 
Fig. A.l: Integral I2/D as function of a and b. 



A.3. Third integral 
Similarly, we can rewrite ^3 as 



7 3 = A 1 ■ h 



o- I dmi I 

a 2 ■ hp r 

v — y — a + 1 J«i 



. l+y-a v-y-a 

d/722 Hlj ' m 2 



dmi m 



v-y-a+l 



-(f) 



(A. 12) 

v-y— a+ 1 \ 



If Eq. A. 10 holds, then from mf » m follows 



h 



A 2 -h 



■ no ■ m f 



v— y— a+ 1 



(v — y — a + 1) • (y - a + 2) 



y-a+2^ 



y-(i'+2 



(A.13) 
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A.4. Deriving the steady-state distribution 

The first and the third integrand I\ and I3 show the same mass 
dependence and can be summed up to 

/1+/3 



A 2 ■ h 

while h is proportional to 
h 



A 2 -h Q 



and therefore 



rrtf 



v-y-a+1 y-a+2 

m f m' , 



h 



h+h 



fflf 



v-a-y+l 



(A.15) 



(A. 16) 



Since the constants of proportionality are factors of order unity 
and m < mf, the integrals I\ + ^3 are much larger than I 2 , there- 
fore, the flux F (to) is proportional to m y ~ a+2 . 

In the case of a steady state, this flux and the downward flux 
of fragments (which is proportional to m 2 ~») have to cancel each 
other. Therefore, the exponents of the mass dependence need to 
cancel out, i.e. 

a = y + £, (A.17) 

which is the slope of the steady-state condition in the fragmen- 
tation dominated regime (Case C in Fig. 1). 

Appendix B: Derivation of the degree of 
homogeneity 

In the following, we will derive the degree of homogeneity (cf. 
Eq. 3) as shown in Table 1. For simplicity, we will drop all 
constant factors and consider only the proportionalities. The 
Brownian motion kernel is a product of the geometrical cross 
71 (a\ + at) 2 and the Brownian motion relative 



section <x 



geo 



velocities (cf. Eq. 32): 



--BM 



oc (a\ + a 2 ) 



m\ + m 2 
ni\ to 2 



m\ ■ (l +6^f (l +6T 1 ) 3 



(B.l) 



=h(9) 



where 6 = mijm\. Thus, the m\ dependence of the kernel gives 
v = \. For 8 <c 1, it follows that h(8) oc Qri and by comparison 
to Eq. 21, we can derive y - — i. 

The relative velocities in the first regime of turbulent relative 
velocities is proportional to \a\ - a 2 \ (cf. Eq. 33), thus, the kernel 
can be written as 



C 1 



OC (fl[ + a 2 ) 

oc mi ■ (1 + ( 



ii) 2 . 



- a 2 \ 



(B.2) 



=/!(«) 



In this case, the function h(6) is different: h(ff) oc 1 for 9« 1, 
and we derive v = 1 and y — 0. 

In the second turbulent regime, the relative velocities are 
proportional to the square root of the Stokes number of the 
larger particle (see Eq. 34). By using a limit representation of 
max(«i, ai), we can write 



v in 



m,,m 2 K («1 +"2) 

oc m\ ■ (l + 



lim 



lim (1 + 03 )2» 

N—>oo 



(B.3) 



Thus, for the limit of small 9, we derive v = | and y = 0. If 
settling is included (see Eq. 31), then v increases by an additional 
factor of 1/6. 



(A.14) Appendix C: Estimating the equilibration time scale 

As shown in the appendix of Birnstiel et al. (2010a), monodis- 
perse growth (i.e. assuming all particles have the same size) pro- 
vides a good estimate of growth time scales. In this approxima- 
tion, the growth rate is given by 



1 



dm 



da 

dt 4np s a 2 dt 



1 



47rp s a 2 T co ii ' 



, (CI) 



where r co u is the collision time scale. Integration of Eq. C.l 
yields the time a particle needs to grow from size a\ to a 2 . For 
Brownian motion and the first two turbulent velocity regimes 
(see Eqs. 32 and 34), we can derive the growth times 



?BM(fll,fl2) 

hu(ai,a 2 ) 
hib(a\,a 2 ) 



1 



2(7T Ps ) 3 



5S d ^k \ 3fim p 
4 i^gjh 



[al-aj] 



7-77- ■ a/ ~t ~ (v«2- v«r) 



2d Qk V 37r \ a 1 



(C.2) 

(C3) 
(C4) 



Here the last two growth times belong to two distinct cases: fn a 
is the time in the second turbulent regime (see Eq. 34) without 
settling effects (i.e., a 2 < a sett ) while fnb includes the fact that 
settling of particles to the mid-plane increases the dust density 
at the mid-plane (denoted by pd) and thus accelerates growth. 

Typically, coagulation starts by Brownian motion growth 
from sub-//m sized particles to sizes where turbulent veloci- 
ties become important (i.e., sizes larger than «bt, see Eq. 37). 
Therefore, we can estimate this time by 



T"BM - ?BM(aO, «BtX 



(C5) 



where ao is the monomer size. 

We will neglect the time needed by particles to grow until 
the second turbulent regime because it is typically much shorter 
than the other time scales involved. If particles in the second 
turbulent regime are already influenced by settling {a\ > a sett ), 
then growth proceeds according to Eq. C.4 and the time is given 
by 

)■ (C6) 

If the size range is entirely below a sett , then the timescale is given 
by 

Tn = flla(fll2.flmax)- (C7) 

For the case that a 12 < a se tt i= a m ax, we need to add both contri- 
butions 

I'll = flIa(fll2,a S ett) + fllb(flsett, flmax)- (C.8) 

By comparison to our time evolving simulations of particle 
growth and fragmentation, we found that 



^"equil — 8 ' (TBM + Tu) 



(C.9) 



estimates within a factor of a few the time at which the distribu- 
tion reaches a steady state. 



